On the influence of time and space correlations on the next earthquake magnitude 
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A crucial point in the debate on feasibility of earthquake prediction is the dependence of an 
earthquake magnitude from past seismicity. Indeed, whilst clustering in time and space is widely 
accepted, much more questionable is the existence of magnitude correlations. The standard approach 
generally assumes that magnitudes are independent and therefore in principle unpredictable. Here 
we show the existence of clustering in magnitude: earthquakes occur with higher probability close 
in time, space and magnitude to previous events. More precisely, the next earthquake tends to 
have a magnitude similar but smaller than the previous one. A dynamical scaling relation between 
magnitude, time and space distances reproduces the complex pattern of magnitude, spatial and 
temporal correlations observed in experimental seismic catalogs. 

PACS numbers: 02.50.Ey,64.60.Ht,89.75.Da,91.30.Dk 
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Since the Omori observation temporal clustering 
is considered a general and distinct feature of seismic 
occurrence. Clustering in space has also been well es- 
tablished and, together with the Omori law and the 
Gutenberg-Richter law [3J] , is the main ingredient of prob- 
abilistic tools for time-dependent seismic hazard evalua- 
tion 0,11,0]. The distribution D(At) of the inter-time At 
elapsed between two successive events is a suitable quan- 
tity to characterize the temporal organization of seis- 
micity. Analogously, the distribution D(Ar) of the dis- 
tance Ar between subsequent epicenters provides useful 
insights in the spatial organization. Both distributions 
have been the subject of much interest in the last years 
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particular, they exhibit universal behavior essentially in- 
dependent of the s pac e region and the magnitude range 
considered d E1HS E3|- Furthermore, the question 
of the existence of correlations between magnitudes of 
subsequent earthquakes has been also recently addressed 
[3 EE IU- In ref.[H EH, Corral has shown that the 
Southern California Catalog exhibits possible magnitude 
correlations that are small but different from zero. How- 
ever, restricting his investigation to earthquakes with At 
greater then 30 minutes, he observes that correlations 
reduce and become smaller than statistical uncertainty. 
Magnitude correlations have been, therefore, interpreted 
as a spurious effect due to short term aftershock incom- 
pleteness (STAI) [2]]]. According to this hypothesis, af- 
tershocks, in particular small events, occurring closely 
after large shocks are not reported in the experimental 
catalog. This interpretation agrees with the standard ap- 
proach that assumes independence of earthquake magni- 
tudes: an earthquake "does not know how large it will 
become" . This has strong implications on the still open 
question of earthquake predictability 22]. On the other 
hand, a recent analysis of the Southern California Cata- 



log [2(| has shown the existence of non-zero magnitude 
correlations, not to be attributed to STAI. These are ob- 
served by means of an averaging procedure that reduces 
statistical fluctuations. A dynamical scaling hypothesis 
relating magnitude to time differences has been proposed 
to explain the observed magnitude correlations. 

In this paper we present a statistical analysis of ex- 
perimental catalogs that confirms the existence of rel- 
evant magnitude correlations. In particular, the anal- 
ysis enlightens the structure of these correlations and 
their relationship with At and Ar. We then intro- 
duce a trigger model based on a dynamical scaling re- 
lation between energy, space and time and show that 
this model reproduces the above experimental find- 
ings. We consider the NCEDC catalog (downloaded at 
http:// /www. 11ccdc.org/11ccdc/. years 1974-2002. South 
(North) lat. 32 (37), West (East) long. -122 (-114)). Sim- 
ilar results are obtained for seismic catalogs of other geo- 
graphic regions. To ensure catalog completeness, we con- 
sider only events with magnitudes m > 3 and take into 
account STAI using the method proposed in ref. [23J . The 
quantities considered are Ati = tj+i— U, Ar,; = |ri+i. — r,;| 
and Ami = "nn+i — tth, i.e. the time, space and magni- 
tude difference between subsequent events. We have also 
evaluated the quantity Am* = m,» — mi where i* 7^ i 
is the random index of an earthquake recorded in the 
catalog. Hence, Am* is the magnitude difference within 
a reshuffled catalog where the magnitude of the subse- 
quent earthquake is independent of previous ones. We 
then consider the conditional probability 
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where N(xo,yo) is the number of couples of subse- 
quent events with both Axi < xo and Ayi < ya 
and N(yo) is the number of couples with Ayi < yo. 



In the following Axi or Ay, will be used to indicate, 
depending on cases, Ar,-, Atj, Ato^ or Am*. Our 
method is schematically presented in Fig.l. Keeping 
m and r fixed, we compute the quantity P(Ato* < 
TOo|Arj < ro) for several independent random realiza- 
tions of the reshuffled catalog, obtaining the distribution 
p [P(Am* < TOo|Arj < ro)]. Taking I0 4 independent re- 
alizations of the magnitude reshuffling, for each given Too 
and ro, we always find that p[P(Am* < too|A^ < ro)] 
is gaussian distributed with mean value Q(roo,ro) and 
standard deviation a(mo,ro). Analogous behaviour is 
obtained for P(Ato* < mo|Ai; < to) and we similarly 
define Q(mo,to) and <t(too,£o)- The relevant quantity is 
5P(m ,yo) = P(Arm < to |A?/; < y Q ) - Q(m ,y ), i.e 
the difference between the value of P(Ato; < mo\Ayi < 
yo) in the real catalog and its mean value in the reshuf- 
fled one. If the absolute value \8P(mo, yo)\ is larger than 
ct(too, yo), significant non-zero correlations between mag- 
nitudes of successive earthquakes exist. In particular, 
a positive value of 5P(mo,yo) > CT (wo,yo) implies that 
the number of couples N(mo, yo) is significantly larger in 
the real catalog with respect to a catalog where magni- 
tudes are uncorrelated. In Fig.l we explicitly compare 
p[P(Am* <m \Ayi <yo)} with P(A mi < m \A yi < 
yo) for mo = and yo = r o = lOfcm or yo = to = lh. 
One clearly observes the existence of non-zero magni- 
tude correlations, since <5P(TOo,ro) ~ 8.3<t(too, To) and 
5P(mo,to) ~ 7.3<t(too, to)- For a deeper understanding 
of the nature of the observed correlations, the above anal- 
ysis has been extended to other values of Too, r and 
to- In Fig. 2 and Fig. 3 we plot the quantities <5P(too,To) 
and 5 P (mo, to) as a function of too for different values 
of ro and to respectively. The error bar of each point 
is the standard deviation a(mo,yo)- We first observe 
that for each value of ro and to and for a wide range 
of too, SP(mo,yo) is strictly positive and significantly 
different from zero. Considering the behavior at fixed 
ro or to, the curve has a peak centered in too < 0, in- 
dicating a crossover from positive to negative correla- 
tions. This can be better enlightened by the derivative 

P'(to ,?/o) = d5P dmo VO%> ~' wri ich represents the proba- 
bility difference for Ato^ = too conditioned to Ayi < yo ■ 
P' '(too, yo) is therefore an estimate of the magnitude cor- 
relation between two subsequent events with Am; = mo. 
Interestingly, for both yo = ro and yo = to (inset of Fig. 2 
and Fig. 3), P'(too,j/o) has the maximum value for too 
in the range [—1,-0.5] and decreases to zero for smaller 
values of too. For too > 0, P'(mo,yo) is always negative 
with the minimum value centered around too G [0, 0.5] 
and going to zero for large too. This implies that, for pos- 
itive too, the probability is larger in a reshuffled catalog 
than in the real one. As a consequence, Fig.s (2,3) clearly 
show that the magnitudes of subsequent earthquakes are 
correlated and, in particular, the next earthquake tends 
to have a magnitude close but smaller than the previous 
one. Furthermore, Fig.s (2,3) indicate that for any fixed 
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FIG. 1: (Color online) The distribution of P(Am* < 
mo\Ari < ro) for ro = WKm and mo = (black curve) 
is compared with P(Am; < 0|Ari < WKm) = 0.543 
(broken black curve). p[P] has a gaussian behaviour 
with mean Q(mo,ro) = 0.473 and standard deviation 
cr(mo,ro) = 0.00825. For ro = WKm and mo = one 
has JP(mo,ro) = 0.07 ~ 8.5cr(mo,ro) strongly support- 
ing the existence of correlations between m, and rrii-i. 
The same conclusion can be obtained by considering 
P(Am* < mo|Ati < to) for to = lh and mo = (red curve). 
It is found P(Am, < 0|At< < lh) = 0.537 (broken red curve) 
whereas p[P] is a gaussian with mean Q(mo,to) = 0.475 and 
standard deviation o"(mo,to) = 0.0085. As a consequence 
SP(mo,t ) ~ 7.3<r(mo,to)- 



to , curves corresponding to different r or t clearly sep- 
arate, showing the existence of correlations between Ato 
and Ar, Ato and At. In particular we observe that the 
larger are ro or to, the smaller are magnitude correlations. 

To better investigate the role of ro and to on magni- 
tude correlations, we consider P(Ar; < ro|ATOi < too) 
andP(Aii < £o|Ato; < too). Following the procedure de- 
scribed for Fig.l, we compute SP(xo,mo) and a(xo,mo) 
for xo = ro,io (Fig. 4). Also in this case, a non-zero 
SP(xo,mo) is the signature of magnitude correlations. 
We observe that, for each value of too, 5P(xo,mo) is a 
decreasing function of ro and to and therefore stronger 
correlations are observed for events that occur closely 
in time and space. More specifically, Fig. 4 shows that 
for smaller values of to and ro, the probability to have 
Ato.; < —2 is about 40% larger in the real than in a 
given reshuffled catalog. 

The above analysis shows that a better description of 
real seismicity can be obtained if correlations between 
time, space and magnitude are properly taken into ac- 
count. As in dynamical critical phenomena where en- 
ergy and time fix a characteristic length scale, similar 
ideas can be used to introduce magnitude correlations 
within standard trigger models for seismicity. In trigger 
models 241 ] . the probability to have the next earthquake 
in the time window [t,t + St], with epicenter in the re- 
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FIG. 2: (Color online) The quantity 8P(mo,ro) as a function 
of mo for ro = 10, 100, 500-Km from top to bottom. For 
each ro and mo the error bar is the standard deviation 
cr(mo,7'o). Data for the Southern California catalog (black) 
are compared with numerical simulations (red). In the inset, 
the quantity P'(mo,ro). 



FIG. 4: (Color online) (a) The quantity <5P(ro,mo) as a 
function of ro for mo = 2, 0.5, from top to bottom, (b) The 
quantity 8P(to, mo) as a function of ro for mo = 2, 0.5, from 
top to bottom. Data for the experimental catalog (black) are 
compared with numerical simulations (red). For each point 
the error bar is the standard deviation cr(a;o,mo). 




FIG. 3: (Color online) The quantity SP(mo,to) as a function 
of mo for to = 1, 10, BOhours from top to bottom. For 
each to and mo the error bar is the standard deviation 
er(mo,io). Data for the Southern California catalog (black) 
are compared with numerical simulations (red). In the inset, 
the quantity P'(mo,io). 



gion [r, r + Sr\ and magnitude in the range [m, m + 8m] 
is given by the superposition 

V(t,r,m) =^P(t - tj,\f - fj\,m,m,j), (2) 



where P(t 



tj, \r 



is the probability condi- 



tioned to the occurrence of an earthquake of magnitude 
rrij, at time tj < t, in the position fj . In the widely ac- 
cepted ETAS model Q, m, rrij, t—tj and \r— fj | are all in- 
dependent quantities and empirical laws are used to char- 
acterize their distributions. Many analytical and numer- 



ical studies show that the ETAS model captures several 
aspects of real seismic occurrence 0, El, 25, 2(| 27, 28j |. 
Nevertheless, because of the assumption of independence 
between m and rtij, 5P{x,y) would be a random fluctu- 
ating function with zero average and standard deviation 
a(x,y), in all cases considered in Fig.s (2-4). Hence, by 
construction, the ETAS model does not take into account 
magnitude correlations and their dependence on time and 
space. 

In order to reproduce the experimental findings, we 
introduce 



ktl0 b( mi -mj) and 



(3) 



which fix two characteristic time scales leading to the 
scaling behaviour with Atij = t i — tj 

P{Atij, \ fi - fj\,rm, raj) = At~f /z H 



(4) 

The exponent 2/z is determined by imposing the condi- 
tion / dfiP^tij, In - fj\,rrii,mj) = B x (s~")> wnere 
the function Pl\{x) must satisfy the normalization con- 
dition J dxH 1 (x) = 1. Following ref. [2(| it is possi- 
ble to show that this normalization removes the problem 
with "ultaviolet" and "infrared" divergences of the ETAS 
model [H. 

In order to simplify the numerical procedure, we con- 
sider a special case of Eq.Q 



H 



At, 



Hi 



i j 



At)( z / VAW \At) /z 



Ho 



(5) 



In the numerical simulation, we generate a synthetic cat- 
alog containing only occurrence times and magnitudes. 



4 



We start with a random event at initial time to = 0, time 
is then increased by one unit and a trial magnitude is ran- 
domly chosen. The i-th event, at time fcj and with mag- 
nitude m t , occurs with a probability X)j<i ^i( T u/^^y ) 
where the sum is over all previous events. In particular 
we use 



Hx{x) = 



A 



e x - 1 



■ 7i 



(6) 



with the parameters k t = 12.7ft., A = 0.21/i _1 , 71 = 0.1 
and b = 0.9, that in ref.[20| are found able to reproduce 
the experimental behavior of magnitude and inter-time 
distributions. To introduce the epicenter location in the 
numerical catalog we use the power law 



H 2 {x) 



B 



x^ + 72 



(7) 



where [i and 72 are fit-parameters and B is fixed by the 
normalization. Other functional forms for H2(x) give 
similar results but with a worse agreement with experi- 

~M for 



mental data. Wc follow the method used in ref . |25|, 
the ETAS model. More precisely, for j — the epicenter 
of the first event in the catalog 7^=0 is randomly fixed in 
a point of a square lattice of size L. For sufficiently large 
lattices, the results are L independent. Next, j is up- 
dated j = j + 1, and the mother of the j — th earthquake 
is chosen among all previous < i < j — 1 events accord- 
ing to the probability Hi (si 2- ) • Once the mother event 
is identified, its epicenter r* and occurrence time t* are 
used to randomly obtain fj from the probability distri- 
bution (tj — t*)~ 2 / z H2 an d assuming space 
isotropy. 

The numerical catalog is analyzed with the same pro- 
cedure applied to experimental data. Numerical results 
for SP(xo,yo) are presented as red circles in Fig.s (2-4) 
for fi = 2.6, k r = 0mh/Km z , z = 3.3 and 72 = 0.1. We 
observe a very good agreement between numerical and 
experimental data. In particular, 5P(nio,yo) (Fig.s(2,3)) 
displays a maximum value localized around the maxi- 
mum of the experimental distribution. Furthermore, also 
the functional form of the decay of SP(x ,m ) for both 

= r o (Fig.4a) and x — t (Fig. 4b) is reproduced. 

The agreement between numerical and experimental 
results, indicates that the scaling relation ((4]) among 
magnitudes, times, and epicenter distances can describe 
the complex pattern of the experimentally observed cor- 
relations. The origin of magnitude correlations, within 
our theoretical approach, has a direct interpretation. Ac- 
cording to Eq.s (|5l6l7p . indeed, at the time t an earth- 
quake of magnitude m with epicenter f has a finite proba- 
bility to be triggered by a previous (rrij , tj , fj ) earthquake 
only if m < mj-(l/b)log((t-tj)/k t ) and 77, < (t— tj) 1 /* 



(Eq. (O). As a consequence, only events occurring 
close in time and space can have a magnitude close and 
smaller, or even larger, that the previous triggering one. 
Magnitude correlations, therefore, become particularly 
relevant within aftershock sequences, when earthquakes 
tend to be very close in time and space. Thus a dynam- 
ical scaling approach, that properly takes into account 
these correlations, can improve existing methods for time 
dependent hazard evaluation. 
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